library(tidyverse)
library(fixest)
library(showtext)
library(vtable)
library(writexl)
library(haven)
library(gridExtra)

setwd('D:\\WorkingPaper-Series\\Township_Light\\Replication_Files\\workdata')

##################################
########## Figure C5 #############
##################################


after_2016 <- read_delim('feature_after_2016.csv')
high_kw <- read_delim('feature_high_kw.csv')

after_2016 <- after_2016 %>% mutate(upper_ci = coef + 1.96*se, lower_ci = coef - 1.96*se)
high_kw <- high_kw %>% mutate(upper_ci = coef + 1.96*se, lower_ci = coef - 1.96*se)

p1 <- after_2016 %>% mutate(after_2016=as.factor(after_2016)) %>% 
  ggplot(aes(type,coef,shape = after_2016, color = after_2016))+
  geom_hline(yintercept = 0, color = 'grey50',linetype =2)+
  geom_errorbar(aes(ymin = lower_ci, ymax = upper_ci),
                width = 0, lwd = .5,position=position_dodge(0.5))+
  geom_point(aes(fill = after_2016),position=position_dodge(0.5))+
  theme_bw()+
  theme(legend.position = c(0.13,0.85),
        legend.background = element_rect(color = 'black', size = 0.1,linewidth = 0.1, fill = 'white'),
        legend.margin = margin(t = -0.4, r = 0.1, b = 0.1, l = 0.1, unit = "cm"),
        text = element_text(family = 'serif',size = 9),
        axis.text =  element_text(family = 'serif',size = 9))+
  labs(y = 'Coefficient', x = '')+
  scale_x_discrete(labels = c('county_dist'='Dist to County GOV','light10'='Initial Light','log_pop'='Initial Population',
                              'log_road'='Initial Road','log_topo'='Topography'))+
  scale_color_manual(values = c('0' = '#279EFF', '1' = '#D80032'), 
                     labels = c('Before 2016','After 2016')) +
  scale_fill_manual(values = c('0' = '#279EFF', '1' = '#D80032'),
                    labels = c('Before 2016','After 2016')) +
  scale_shape_manual(values = c('0' = 21, '1' = 24), 
                     labels = c('Before 2016','After 2016'))+
  guides(color = guide_legend(title = ""),shape = guide_legend(title = ""),fill = guide_legend(title = ""))

p2 <- high_kw %>% mutate(high_kw=as.factor(high_kw)) %>%
  ggplot(aes(type,coef,shape = high_kw, color = high_kw))+
  geom_hline(yintercept = 0, color = 'grey50',linetype =2)+
  geom_errorbar(aes(ymin = lower_ci, ymax = upper_ci),
                width = 0, lwd = .5,position=position_dodge(0.5))+
  geom_point(aes(fill = high_kw),position=position_dodge(0.5))+
  theme_bw()+
  theme(legend.position = c(0.13,0.85),
        legend.background = element_rect(color = 'black', size = 0.1,linewidth = 0.1, fill = 'white'),
        legend.margin = margin(t = -0.4, r = 0.1, b = 0.1, l = 0.1, unit = "cm"),
        text = element_text(family = 'serif',size = 9),
        axis.text =  element_text(family = 'serif',size = 9))+
  labs(y = 'Coefficient', x = '')+
  scale_x_discrete(labels = c('county_dist'='Dist to County GOV','light10'='Initial Light','log_pop'='Initial Population',
                              'log_road'='Initial Road','log_topo'='Topography'))+
  scale_color_manual(values = c('0' = '#279EFF', '1' = '#D80032'), 
                     labels = c('Low KW Freq','High KW Freq')) +
  scale_fill_manual(values = c('0' = '#279EFF', '1' = '#D80032'),
                    labels = c('Low KW Freq','High KW Freq')) +
  scale_shape_manual(values = c('0' = 21, '1' = 24), 
                     labels = c('Low KW Freq','High KW Freq'))+
  guides(color = guide_legend(title = ""),shape = guide_legend(title = ""),fill = guide_legend(title = ""))